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Abstract 

Motivated by recent developments in ion trap design and fabrication, we investigate the stability of ion 
motion in asymmetrical, planar versions of the classic Paul trap. The equations of motion of an ion in such a 
trap are generally coupled due to a nonzero relative angle between the principal axes of RF and DC fields, 
invalidating the assumptions behind the standard stability analysis for symmetric Paul traps. We obtain 
stability diagrams for the coupled system for various values of d, generalizing the standard q-a stability 
diagrams. We use multi-scale perturbation theory to obtain approximate formulas for the boundaries of the 
primary stability region and obtain some of the stability boundaries independently by using the method of 
infinite determinants. We cross-check the consistency of the results of these methods. 

Our results show that while the primary stability region is quite robust to changes in 9, a secondary 
stability region is highly variable, joining the primary stability region at the special case of 6^ = 45°, which 
results in a significantly enlarged stability region for this particular angle. 

We conclude that while the stability diagrams for classical, symmetric Paul traps are not entirely accurate 
for asymmetric surface traps (or for other types of traps with a relative angle between the RF and DC 
axes), they are "safe" in the sense that operating conditions deemed stable according to standard stability 
plots are in fact stable for asymmetric traps, as well. By ignoring the coupling in the equations, one only 
underestimates the size of the primary stability region. 



1 Introduction 

In the quest to miniaturize ion traps for quantum information processing (QIP), many of the recent designs 
being explored are planar versions of linear Paul traps. [l}|3] Since the electrodes all lie in a single plane, 
these traps can be constructed using VLSI microfabrication, which offers great scalability and potential to be 
integrated with other useful on-chip components such as mirrors, fiber ferrules, and cavities. The ions trapped 
by a surface trap are cooled by laser beams that are commonly aligned parallel to the trap surface. In order 
to cool the ions efficiently, the three principal axes of the trap pseudopotential must have nonzero components 
along the laser direction [i]. In particular, none of the principal axes can be vertical, if the cooling beam is 
parallel to the horizontal trap surface. The necessary tilt in the principal axes relative to the beam direction is 
typically realized by using asymmetric electrode designs and/or setting the voltages in an asymmetrical manner. 
In general, such asymmetry, in addition to introducing the desired tilt, introduces a relative angle between the 
principal axes of the RF and DC fields. (See Figure [l]) 
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Figure 1: The cross-section of an example surface trap and the associated electric fields. Figure (a) shows the 
RF electrodes (in red) and the DC control electrodes (in yellow) on a substrate and the total pseudopotential. 
The DC electrodes are used to trap the ion axially. The trap axis is orthogonal to the page, and the trap center 
(the location of the ion) is at the minimum of the pseudopotential, shown in blue. Figures (b) and (c) show 
the RF and DC potentials near the trap center, as indicated. The potentials were calculated numerically using 
the trap geometry. As can be seen, the RF and DC principal axes are not aligned perfectly, and have a small 
angle between them. Due to this nonzero relative angle, the ion motion is coupled in the two radial directions, 
independent of the orientation of the coordinate axes. 

When the angle between the RF and DC principal axes is nonzero, the classical equations of motion of an 
ion near the trap center are given by a coupled version of the Mathieu equation. The stability properties of such 
a coupled system cannot be obtained from the classical stability analysis of Paul traps, which assumes that the 
equations are decoupled. Thus, it is not clear, a priori, that the operating conditions obtained by resorting to 
the stability analysis of symmetric Paul traps will result in stable ion motion in asymmetric surface traps, as 
well. One needs to do the stability analysis for the asymmetric case from scratch, and obtain the corresponding 
stable operating conditions. 

In this paper, we generalize the standard q-a stability diagrams for symmetric Paul traps to obtain stability 
diagrams for asymmetric surface traps, or more generally, for trap designs and operating conditions that result 
in a relative angle between the radial axes of RF and DC fields. We also obtain approximate formulas for 
the boundaries of the primary stability region, generalizing the formulas for the symmetric, decoupled case. 
These results give the stable operating conditions for asymmetric surface traps, and can serve as a reference for 
experimentalists trying to select regions of stability for ion motion, instead of having to rely on the untested 
assumption that the diagrams for symmetric Paul traps are still applicable. 

We approach the problem by both numerical and analytical techniques. On the numerical side, for a given 
set of parameter values describing the operating conditions of the system (such as the generalizations of the 
q and a parameters of the symmetric Paul trap), we obtain a basis set of numerical solutions, and by using 
results from Floquet theory, decide whether the system is stable or unstable under the given conditions. By 
scanning a range of values of the parameters and collecting stability data over a region in the parameter space, 
we obtain the relevant stability diagrams. Our analysis confirms the stability theory prediction [5j[6] that as 
the parameters are varied, instabilities develop when the eigenvalues of a certain solution matrix collide on the 
unit circle. As an independent check, we utilize the method of infinite determinants [7 to directly obtain some 
boundaries of the stable regions, and show that these curves agree with the results from Floquet theory. 

On the analytic side, we utilize multi-scale perturbation theory to obtain approximations to the curves 
bounding the primary stability region [8j|9]. These approximate results take the form of formulas that relate 
the parameters describing the system, such as the generalized q and a parameters mentioned above, and the 
angle between the RF and DC principal axes. Certain novelties of the coupled multi-variable case complicate 
the analysis, and we use a hint from our numerical results to pick the relevant curves. 

We show that our analytical results well approximate the boundaries of the stability domains in their regimes 
of applicability, except for the special case of a 45 degree tilt between the RF and DC principal axes. We comment 
on this special case, for which the fundamental stability domain is significantly enlarged. 

The paper is organized as follows. In Section [2] we set up the general equations of motion of the coupled 
system, restricting our attention to the two- variable case. In Section |3j we give a discussion of the stability of 
periodic Hamiltonian systems, and then apply the formalism to the specific case of the coupled Mathieu system. 
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obtaining stability diagrams via numerical solutions of the equations. In Section |4j we discuss an alternative 
method, the method of infinite determinants, which is capable of obtaining some of the stability boundaries 
directly, and demonstrate its consistency with the results of Section [3] In Section [5] we turn to the method of 
multiple scales, and after a general discussion, apply this method to the coupled Mathieu system relevant for 
surface traps. This analysis results in approximate formulas for the boundaries of the primary stability region, 
which we check against the numerical results of Section |3] Finally, in Section [6] we summarize our findings, and 
discuss the practical implications for ion trap design and operation. 



2 Equations of motion. 

Let us begin by writing the general equations of motion of an ion near the center of an asymmetrical trap. We 
will assume that the oscillating (RF) and static (DC) electric fields in the ion trap have a coincident zero at 
the trap center, which we take to be the origin of our coordinate sysetm, x = {x, y, z) = (0, 0, 0). We will work 
in the harmonic approximation and treat potentials as second order in the displacements from the origin, and 
forces (or electric fields) as first order. Since the electric fields vanish at the origin, the first order terms in 
the expansions of the potentials vanish. Denoting the potential energy of an ion hy U = eV, where e is the ion 
charge, and choosing the zero of U so that C/(0) — 0, we have, up to second order in the displacements, 

U{x,y,z,t) = U^'^ {x,y,z) cos Lot + U^'^{x,y,z) (1) 
= ^^x,a;jC/«Fcos(u;i) + ^^a;,x,C^DC^ (2) 

where Xi,X2,X3 stand for x,y,z, respectively, lu is the RF angular frequency, and the Ufj^ and U^'~^ are the 
(symmetric) matrices of second derivatives of the RF and DC potential energies. 
The equations of motion are given as, 

(fxi dU 

j j 

Defining a new time variable t hy ujt ~ 2t and denoting the derivatives with respect to t by dots, we get. 



Xi + AijXj + 2 QijXj cos 2t — , (5) 

3 3 

where. 



are the multi-variable versions of stability parameters a and q of the Mathieu equation 10 . stiffness matrix. 
Qij are traceless, 

^A,, = 0, J2q,,=0. (7) 

i i 

The decoupled case. One can diagonalize Uj^'^' by an orthogonal transformation S, 

Xi = SijXj, (8) 
j 

and obtain, 

U^^ii, yrz)-\ + + Uf^ii) , (9) 



^We work in the quasi-static approximation, where the RF field is described by an instantaneous electrostatic potential. 



3 



where the diagonal matrix Uij is given by Uij = '^i^i SikSjiUki- If the RF and DC principal axes coincide, 
the same transformation ^ also diagonalizes , and we get (dropping the tildes, in order to simplify the 
notation) , 

U= '^{U^^xl + U^^xl + UWxl)cosicot) + 



DC ^2 



22 



u. 



DC^2 



33 



0- 



This is the classical case of a symmetric Paul trap. The equations of motion resulting from this potential are 
decoupled; Equations ([s]) and (|6| become. 



2qi cos 2T)xi = , 



(10) 



where. 



(11) 



Equation ( 10 ) is known as the single- variable Mathieu equation, and by using the results of the classical Mathieu 



stability analysis on each component separately, one can obtain the regions of joint stability. This gives the 



standard a-q stability plots for the Paul trap 11 12 . 



Reduction to two dimensions. In most surface trap designs, the RF electrodes are long, and the total RF 
field has non-vanishing components only in the radial directions, which we denote by x and y. If, in addition, 
the electrodes and voltages are symmetric around the z = plane, which is commonly the case, the axial (z) 
motion of the ion is decoupled from the x and y motions, and is simple harmonic (since the axial force in this 
case is given by a DC field linear in the displacement z). The x-y equations of motion arc still of the form ([5|, 

Xi + ^ AijXj + 2 ^ QijXj cos 2t = , (12) 

j J 

but the matrices A and Q are now 2x2. When the DC field in the z direction is confining. Gauss's law implies 
that the DC field has an anti-confining radial component. In other words, the trace of U^^'^ , as restricted to the 
x-y plane, must be negative. Due to (|6| the same is true for the matrix A. Similarly, since the z-component of 
the RF field is assumed to be zero identically, Gauss's law enforces the 2x2 version of the matrix to be 
traceless. Hence, Q is also traceless. 



A convenient coordinate system and parametrization. Given the general equations (12 1, one can in- 
vestigate the stability properties of the system in terms of the entries of the (now 2x2) matrices Q and A. 
However, in order to obtain stability plots that reduce, in the decoupled limit, to the familiar q-a stability plots 
of symmetric Paul traps, we will fix the relative magnitudes of the entries of Q and A, and vary their overall 
scales. In terms of the actual operating conditions of a trap, this amounts to fixing the trap geometry and the 
ratios of the DC electrode voltages, and changing the overall scale of the voltages on DC and RF electrodes. 
The matrices A and Q are symmetric since they are related to the symmetric matrices Uf^^ and Uj^'^ of 
through ([6|. Thus, we can diagonalize at least one of A or Q by a suitable choice of our coordinate axes 
X and y. Let us assume that A is diagonalized, and write it in the form, 

^ = 4l (13) 



X -a ^ 

where a and a are constants to be determined by the electrode geometry and the DC voltages. As argued above, 
A must have negative trace due to Gauss's law and the fact that the ion is confined along the z-axis. Thus, we 
must either have a > and a > 1, or a < and a < 10 

Recall that U^^ is symmetric and traceless. If we were to use a coordinate system (i, y) in which U^^ is 
diagonal, the amplitude of the RF potential energy would have the form, 

U^''ii,y) = lurii'-f), (14) 



^Having a < and a < would correspond to the DC potential being anti-confining along both radial axes. Although this is 
possible, it is not frequently encountered in the designs used in practice. 
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and Qij, being related to Uf^ through (|6l would be given as, 



Q 



1 

-1 



(15) 



where g is a free parameter. The two coordinate systems, (i, y) (in which Q is diagonal) and {x, y) (in which 
A is diagonal) are related by a rotation. 



X 

y 



X cos tl + y sm o 
—xshi6 + j/cos( 



(16) 
(17) 



where 9 is the angle of rotation. Substituting these in (14 1, we get the RF potential energy in terms of the 
coordinates along the DC principal axes. This gives. 



[/R-F (x, y) = i U^'^ (a;2 cos 20 - cos 26 + 2xy sin 29) cos 2t . 



Using ^ and we get the resulting Q matrix as. 



cos 20 sin 29 
sin 26* -cos 26* 



Finally, using (13) and (19), the equations of motion (12) become 



X + ax + 2q{cx + sy) cos 2t = 
y — aay + 2q{sx — cy) cos 2t = 0, 



(18) 



(19) 



(20) 
(21) 



where for brevity we replaced cos 29 and sin 20 by c and s, respectively. The classic case of symmetric Paul traps 
where 9 — is recovered by setting s = and c — I. Below, we will get q-a stability plots for various values of 
a and 9. As mentioned above, fixing a and 9 and varying q and a correspond to fixing the trap geometry and 
the ratios of the DC voltages, and changing the overall scale of the RF (g) and DC (a) voltages. 



We next turn to the investigation of the stability properties of the coupled equations (20)-(21) 



3 Stability of the coupled Mathieu system 
3.1 Periodic systems and stability 

We begin by reviewing some general aspects of linear, periodic systems and their stability. The discussion is 



perhaps a bit abstract, but the bottom line is the following: The stability properties of the system (12) is 



determined by first obtaining a basis set of solutions over one period of the RF field, and then inspecting the 
eigenvalues of the matrix formed by joining these fundamental solutions. 



Equivalent first order system. Let us first rewrite ( 12 ) as an equivalent, first order system by defining the 
velocity components as new variables. Letting, 



x 
x 



(22) 



where x — {x,y) is the radial position vector, we can rewrite the radial equations of motion (12) as, 

u = G(t)u, 

where, 

G(r) = 



(23) 
(24) 



I 

-2Qcos2t-A 

I denoting the 2x2 identity matrix. The matrix G is periodic in time with period T = tt. 
Although our primary interest will be in the two-dimensional case 

for which u is a 4-dimensional vector, much of what we will say below will be valid for a general Hamiltonian 



system (23) with a periodic G(r). 
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Fundamental solution matrix. In order to explore the long-time stability properties of ( 23 ) , we first obtain 



a set of fundamental solutions Ui(T) that form a basis for the space of all solutions. We take the initial value 
Uj(0) of the ith fundamental solution to be the ith column of the 4 x 4-dimcnsional identity matrix. In other 
words, we take the zth component of Ui(0) to be 1, all the other components to be 0. Since the system (23) is 
linear, a solution with arbitrary initial condition u(0) = Uq can be obtained as an appropriate linear combination 
of these fundamental solutions, 



u(r) 



E 



We combine the column vectors Ui{T) into a 4 x 4 "fundamental solution matrix" 

U(r) = [ui(r)u2(T) U3(t) U4(r)] . 



Since each column of the U satisfies (23), U itself satisfies, 



with the initial condition, 



The general solution ( 25 ) can be written as 



U = G(t)U, 
U(0)=/. 

u(r) - U(r)uo . 



(25) 

(26) 

(27) 
(28) 

(29) 



If the system (23) were autonomous, that is, if G(t) were independent of time, it would be possible to treat the 
fundamental solution matrix U(t) as a "time translation operator". An analogue of (29) would be valid for an 
arbitrary initial time tq — given initial conditions u(ro) = Uq, one could obtain the solution at time tq + r by 
applying the map U(t). However, since the system under consideration is not autonomous, this is not the case. 
The matrix U only gives time translations from the moment r = 0, as in (29). 



Mapping at a period. Although the system (23) is not invariant under arbitrary time translations, it is 



invariant under translations by one period T = tt of the periodic matrix G(r). Using this fact, it is possible to 
prove that, 



U(toT) = [U(T)]' 



(30) 



Thus, repeated applications of the matrix U(T), called the matrix for the "mapping at a period", give relevant 
information about the long time behavior of the system ( 23 ). In particular, it can be shown that the equilibrium 
solution of ((23|, u(r) = 0, is stable if and only if the zero vector Uo = 0, is stable under successive applications 
of U(T) ^|6]"F^ If there is a vector v such that [U(r)]'" v for m = 1, 2, ... is an unbounded sequence of vectors. 



then there will be an initial condition of the system ( 23 ) arbitrarily close to origin, for which the solution will 
grow unboundedly. 



This connection between the stability of (23) and the stability under successive applications of the matrix 
U(r) opens the door to testing for the stability of (23) by investigating the eigenvalue-eigenvectoi]^ decompo- 
sition of the matrix U(T) [5][6]. In particular, the equilibrium is unstable if U(T) has an eigenvalue A with 
|A| > 1. If the eigenvalues of G(T) are not repeated, and if all have magnitudes less than or equal to 1, then 
the equilibrium is stable. 



Eigenvalue spectrum of the mapping at a period for Hamiltonian systems. While the results men- 



tioned above are valid for any periodic G(r), the specific case of (24) is obtained from a Hamiltonian system. 



and the theory of classical mechanics imposes certain conditions on the eigenvalues of the relevant U(r). Time 
evolution in Hamiltonian systems is a canonical transformation, and this forces U(T) to satisfy |5], 



u^(r)JU(r) = J, 



(31) 



^ In the sense that small deviations from the vector remain small under successive applications of U(T). 
*or more generally, Jordan 
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where J is an antisymmetric matrix given in terms of the 2x2 identity matrix I as, 



I 

1 



(32) 



Using (31 ), it is possible to show that if A is an eigenvalue of U(r), then, so is A ^. Since the coefHcients of the 



characteristic polynomial of U(T) are all real, we see that A, the complex conjugate of A, is also an eigenvalue. 
Thus, the eigenvalues of U(T) come in four kinds of groups: 

• 4-tuples A, A, 1/A, 1/A, where A 7^ A, |A| ^ 1. 

• Real pairs, A, 1/A, where A = A, |A| 7^ 1. 

• Pairs on the unit circle. A, A, where 1/A = A 

• Lone eigenvalues, A = 1. 

The multiplicities of a given 4-tuple (or a given pair in the real or unit length cases) are the same. 

It follows from this list of possibilities that the trivial solution u = is stable under successive applications 
of U(r) only if all the eigenvalues of U(T) are on the unit circle; the remanining cases must have at least one 
eigenvalue with magnitude larger than 1. If, in addition to being of unit magnitude, all the eigenvalues are 
distinct, then the system is necessar ily stablejf] 

As we change the operating conditions of the trap, the parameters describing the system (23), i.e., the 
entries of the matrices Q and A change. This results in a change in the solutions, and hence, the eigenvalues 
of the mapping at a period, U(r). Suppose we start with a set of distinct eigenvalues on the unit circle, 
corresponding to a stable operating condition of the trap. The above list of possibilities ensures that as we 
tweak the parameters (the voltages), the eigenvalues of U(T) can move off the unit circle only if they "collide" 
on the unit circle first, meaning that instabilities develop through collisions of eigenvalues of unit magnitude. 
We will verify this prediction below, for the case of the coupled Mathieu equations]^ 



3.2 Application to the Mathieu system 



We next apply the general results reviewed above to the coupled Mathieu system, (20l-(21). We proceed as 



follows. For given, specific values for a and 0, we loop over pre-chosen ranges of q and a. For each pair of values 



of q and a, we solve the equations (20)-(21) numerically, for four different initial conditions; we set one of a;(0), 
2/(0), i(0), y(0) to one, and all the others to zero. After obtaining the four solutions over one period of the 
oscillating RF field (i.e., in terms of the time parameter r, over a time period of T = tt), we obtain U(T), the 
"mapping at a period", by joining the four column vectors u.^. 



i^(T) 



(33) 



We then obtain the eigenvalues of this square matrix, and determine the stability properties of the system by 
checking the magnitudes of the eigenvalues. According to the discussion above, if all eigenvalues are distinct 
and are of unit magnitude, then the system is stable, and if there is an eigenvalue of magnitude larger than one, 
the system is unstable. We go beyond distinguishing between instability and stability, and also note whether 
the system is partially stable, having one pair of eigenvalues on the unit circle, and one pair away from it. While 
partial instability means instability in real life (since exciting only the stable subspace of solutions is impossible 
due to noise) , the distinction between partial instability and full instability will help us below in our discussion 
of the multi-scale perturbation analysis. 

Looping over values of q and a and recording the stability properties of the system for each pair of values, 
we obtain a stability diagram for given values of a and Q. We present two such plots below in Figures [2] and [5] 



and 



^If some of the eigenvalues are degenerate, one has to consider the Jordan decomposition [5j|6j. 

®Note that not all such collisions of eigenvalues result in instabilities. A more detailed discussion of these issues is given in [H] 
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the first one being an example of the classical, decoupled Mathieu system, and the second one being an example 
of the coupled case. As mentioned above, we go beyond the binary classification of stable vs. unstable points, 
and indicate regions of partial stability by using a shade of gray. 

In Figures |3]|4] and |6][7) we show how the eigenvalues evolve as one moves on the stability plot, moving across 
boundaries of stability. As predicted in our discussion above, changes in the degree of stability are accompanied 
by collisions of the eigenvalues on the unit circle. We will give a more comprehensive set of stability plots 
after we present an approximate analytical method (the method of multiple scales) for obtaining the stability 
boundaries. 

0. 0.5 I. 1.5 2. 

2.16 



1.62 



1.08 



0.54 



0. 



-0.54 



-1.08 



-1.62 



-2.16 

0. 05 1. 1.5 2. 

Figure 2: Stability plot for the decoupled system; ~ 0° , a = 1/2. The stability for each value of q (a;-axis) 
and a (y-axis) is indicated by the gray level, which represents the number of unit length eigenvalues of U(T), 
the "mapping at a period" . Black indicates complete stability, with all four eigenvalues being on the unit circle. 
Gray indicates partial stability, with 2 eigenvalues on the unit circle, two away from it. White indicates full 
instability, with all four eigenvalues being off the unit circle. The two vertical lines will be referred to in Figures 
Handll 
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Figure 3: The evolution of the eigenvalues of the "mapping at a period" , as one moves along the first vertical line 
(near q = 0.5) in Figure [5] Each row of four plots corresponds to a change in the number of stable eigenvalues 
as one moves along the first vertical line in Figure [2] As described in the text, a change in the number of stable 
eigenvectors is effected by a "collision" of eigenvalues on the unit circle, with eigenvalues on the unit circle 
leaving, or eigenvalues off the unit circle getting on the unit circle. This behavior is confirmed in these figures. 
In some plots, one of the four eigenvalues is outside the region shown, so only three eigenvalues are seen. Note 
that in this decoupled case, all collisions happen on the real line. 
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Figure 4: Similar to Figure [3) the evolution of the eigenvalues of the "mapping at a period", as one moves along 
the second vertical line (near q — 1.5) in Figure [2[ Once again, for this decoupled case, collisions happen only 
on the real line. 
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0. 0.5 1. 1.5 2. 




0. 0.5 1. 1.5 2. 



Figure 5: Stability plot for — 6.4°, a — 1/2. Due to the nonzero angle between RF and DC principal axes, 
the X and y motions are now coupled, and this case cannot be investigated by the standard, single-variable 
Mathieu techniques. As in Figure [2] the number of stable eigenvalues of U(T) is indicated by the gray level, 
black corresponding to complete stability and white corresponding to complete instability. 
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Figure 6: Similar to Figure [3j but for the coupled system whose stabihty plot is given in Figure [Sj Each row 
shows the evolution of the eigenvalues of the "mapping at a period" as one moves along the first vertical line 
(near q — 0.5) in Figure [sj near points where the number of stable eigenvalues changes. Such changes are 
represented by changes in the darkness of the q-a area plot in Figure [5] The transitions on this first vertical 
line (near q = 0.5) all have analogues in the decoupled case, and the "collisions" of eigenvalues still happen on 
the real line only. 
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Figure 7: Similar to Figure |4] but for the coupled system of Figure [5] Each row shows the evolution of the 
eigenvalues of the "mapping at a period" as one moves along the second vertical line (near q = 1.5) in Figure 
[5j near points where the number of stable eigenvalues changes. This time, we encounter a so-called "combined 
resonance^^ near a — 0.9, which doesn't have an analogue in the decoupled case. This is represented as a collision 
of eigenvalues on the unit circle, away from the real line. 



4 The infinite determinant method 



The "infinite determinant" approach to the stability analysis of the Mathieu equation 13 consists of substituting 
a modified Fourier expansion (a Floquet expansion) into the equations of motion and obtaining an infinite 
set of linear equations for the Fourier coefficients. These equations have nontrivial solutions only when the 
determinant of a certain infinite rank matrix (involving the parameters q, a, a, 9) vanishes. In practice, the 
infinite determinant is replaced by the determinant of a large-rank, finite matrix, and setting this determinant 
to zero gives approximate results. We will describe this method in the setting of the single variable Mathieu 
equation, explain how it can be generalized to the multi-variable case, and present plots for certain "simple" 
stability boundaries of the coupled system that can be obtained from this method in a straightforward way. 
These boundaries are in agreement with the results of the numerical approach described in the previous section. 
Using Floquet 's theorem, we look for a solution to, 

x+ (a + 2(7COs2r)a; = 0, (34) 

of the form, 

oo 

x{T)^e'-^ &ne'2"^ (35) 



13 



Here, the infinite sum represents a periodic function with period equal to the period of the sinusoidal term in 



(34), and the exponential term in front determines the long time stability of the solution. The system is stable 
when V is purely real, and unstable when v has a nonvanishing imaginary partQ Substituting (35) into (34) 



and shifting the summation index in two of the terms, we get, 

DO 

^^.r J2 [(-(;. + 2n)2 + a)6„ + g(6„„i+6„+i)]e'2„r^Q^ (3g) 

n— — oo 

Setting the coefficient of each basis functions e*^""^ to zero, we obtain an infinite set of equations for the 
coefficients &„. This set of equations has a solution only if the determinant of the infinite matrix of coefficients 
vanishes, i.e., detB = 0, where, 

: m ~ n 

: m = n ±1 (37) 
: Otherwise. 




The determinant of the matrix (37) does not converge as it standsj^ however, in order to extract numerical 
results, we will work with a finite-size version of the matrix, and the convergence issue will not be of practical 
concern for the sizes and the numerical precision we will be working with. 

Setting the determinant of B to zero gives an equation that relates q, a, and the exponent v. Once again, for 
given q and a, we can extract the growth factor, this time by solving the resulting algebraic equation. Checking 
whether the equation for has an imaginary solution allows us to deduce the stability/instability of the system 
for the given values of q and a. Looping over values of q and a and repeating the analysis would give the 
relevant stability plots. For the single variable case (or for the decoupled, two- variable case), it is possible to 
extract the stability boundaries directly, without looping over q and a and finding the transitions from stability 
to instability — it turns out [t] that for this case, setting the growth factor e'"'^ to ±1 gives equations that relate 
q and a on the stability boundaries. 

The generalization of the infinite determinant approach to the multi-variable case is straightforward, and 
involves replacing various scalars by vectors/matrices. One still gets a determinant equation relating a, q and 
but in this case, not all stability boundaries can be obtained directly, without looping over q and a. Setting 
e^'^'^ to ±1 as in the single variable case gives the stability boundaries that are related to the so-called "natural 
resonances" , but there is another set of stability boundaries, namely, those related to "combined resonances" , 
which cannot be obtained by this method. In 7 , a pragmatic but unrigorous approach is proposed for obtaining 
the boundaries corresponding to the combined resonances, but we do not follow this procedure here. Instead, 
we only obtain the boundaries that are related to natural resonances, and compare the results with those we 
get from the numerical approach described in the previous section. 

In Figures [slim we give the regions of complete stability obtained by the methods of the previous section, and 
the stability boundaries due to the "natural resonances" , obtained by the infinite determinant method described 
here. We choose a = 1/2, and show plots for 4 different values of 6. As can be seen, the boundaries obtained 
by the infinite determinant method bound the primary stability region accurately, except for the special case of 
6 = 45° Smaller secondary stability regions, which move as the angle 9 changes, also have two boundaries that 
are given accurately by the infinite determinant analysis. However, these regions also have a third boundary 
that is invisible to the method employed here, which is capable of getting only the boundaries due to natural 
resonances. 

'^Positive imaginary parts result in decaying solutions, however, our discussion in Section [s] implies that such solutions are 
accompanied with solutions which have u values with negative imaginary parts, i.e., t hose that grow unboundedly in time. 

**We could remedy this by obtaining an alternative set of equations equivalent to \3fi\ by dividing each equation derived from 
l |36| by the corresponding diagonal entry, Bnn- Working with the matrix resulting from such an alternative set of equations results 
in a convergent determinant. 

^We will give a more comprehensive set of stability plots in the following sections, where the behavior around 9 = 45° will 
become more transparent. 



14 



0.0 0.5 1.0 1.5 2.0 

Figure 8: Stability boundaries obtained from the infinite determinant method (the curves) shown together with 
the primary stabihty region obtained by numerical analysis (dark region), for the decoupled, 9 — 0° case, and 
a — 1/2. The curves obtained from the infinite determinant method indeed form boundaries for the primary 
stability region. 



1.0 - 



0.5 -I 



0.0 



-0.5 - 



-1.0 -I 
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Figure 9: Stability boundaries obtained from the infinite determinant method (curves) shown together with the 
stability regions obtained by numerical analysis (dark regions), for 9 — 12°, a = 1/2. The curves give all the 
boundaries of the primary stability region, but not all the boundaries of a small, secondary stability region. 



15 
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0.0 0.5 1.0 1.5 2.0 

Figure 10: Similar to Figures [8] and |9j for 9 = 32°. As in Figure [9j the curves obtained from the infinite 
determinant method form boundaries for the primary stability region, but do not give all the boundaries of two 
small, secondary stability regions. 
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0.0 0.5 1.0 1.5 2.0 

Figure 11: Similar to Figures [8| and [9l \lO\ for the special angle 6 — 45°. The red curves obtained from the 
infinite determinant method fail to give the boundaries of the enhanced stability region, which is formed by the 
joining of the primary stability region with the secondary regions. See the following sections for more on the 
behavior around this special angle. 

5 Method of multiple scales 

For a given set of parameters, the Floquet analysis of Section |3.2| gives the stability regions of the coupled 
equations of motion numerically. While such results are useful, for asymmetric trap design and characterization, 
it would be beneficial to have at least approximate analytical formulae for the stability boundaries. Such 
approximate formulae would enable trap designers to quickly check the stability regions for a particular trap 
geometry, and would guide the trap design. Analogous results for the case of symmetric traps are well-known 
and useful. 

One possible approach to this problem is to use straightforward perturbation theory of ordinary differential 
equations, which involves treating a certain parameter (g, in our case) as small, and expanding the unknown 
solution into a power series in terms of this parameter. This gives separate equations for each power of the 
perturbation parameter, which are then solved order by order. While useful in many settings, this approach fails 
for a variety of problems due to the appearance of so-called "secular", or "resonant" terms, which grow in time 
in an unbounded manner even when the exact solution is known to be bounded. This behavior invalidates the 
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assumptions in the perturbation analysis after a certain time period. A straightforward perturbation expansion 
of the Mathieu equation results in such behavior. 

Multi-scale perturbation theory, or the method of multiple scales, is a more sophisticated technique capable of 
avoiding such secular terms by expanding not only the solution, but also the independent parameter (typically, 
time) to a series in the small parameter. This results in a set of partial differential equations (PDEs) that 
collectively represent the ordinary diffential equation under consideration. The solutions of these PDEs are 
then restricted by imposing the condition that secular terms do not appear. The absence of secular terms allows 
the method of multiple scales to give uniformly accurate results for cases where standard perturbation theory 
gives results valid only for short time intervals. In this Section, we will use this technique to obtain approximate 



stability boundaries for the coupled system, (20l-(21). Since the calculations are a bit tedious, for readers only 



interested in the final results, we note that the approximate formulas for the stability boundaries are given in 



Equations (104), (105), (132) and (137) 



The method of multiple scales has been applied to a wide variety of problems in physics, engineering, 
and applied mathematics. There are various conventions followed in the literature; below, we will use the 
approach described in [o] and [s]. For an elementary introduction to the method of multiple scales, see [9l[l4]. 
For applications to the coupled Mathieu system, see (8[[T5]. Note that the analyses in [s] and [l5] are not 
immediately applicable to our case, since, as opposed to the cases covered in these references, some of the 
natural frequencies for the Paul traps we are considering are imaginary (i.e., one of the transverse directions is 
anti-confining) . 



5.1 Single- variable Mathieu equation 

Let us first demonstrate the use of multi-scale perturbation theory on the relatively simple case of the single- 
variable Mathieu equation. After this example, we will work out the multi-scale perturbation analysis of the 



coupled system (20)-(21) describing asymmetric surface traps. 



Preliminary discussion. We would like to obtain approximate formulas for a{q) on the stability boundaries 
of the Mathieu equation, 

i-t- (a + 2gcos2T)a; = 0. (38) 

This equation has multiple regions of stability separated by regions of instability; we will focus on the primary 
stability region near the origin of the a-q plane. We will treat g as a small variable, and perform an expansion 
in its powers. 

We begin the multi-scale perturbation analysis by introducing a set of slowly-varying time variables (the 
"multiple scales"), 

To = T, Ti = qr, T2 = g^r, .... (39) 



We replace the dependent variable x with a function of the T^s, 

x{t) = x{To,Ti,...). 

and expand it into a series in g, 

x{To,Ti, . . .) = a;o(7o,Ti, . . .) -I- qxi{To,Ti, . . .] 



(40) 



(41) 



The chain rule gives the derivative with respect to r in terms of partial derivatives with respect to the time 
variables T,-, 



d 



d 



d 
dfi 



2 d 



(42) 



Substituting (40), (41) and (42) into (38 1, one gets the equations for multi-scale perturbation theory. These 



equations give approximate solutions to (38) that are valid for small values of q, unless a is close to certain 



critical values. These critical values are of special interest to us, since they are the locations where the stability 
boundaries intersect the q = axis. Various approaches exist for dealing with these critical points; we will 
follow the technique used in where one expands a into a series, as well, around a critical point oq, 



ao+aiq + a2q + 



(43) 
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In this approach, one obtains the stability boundaries by first getting conditions on oi, 02,... that ensure 
that the approximate solutions are bounded. These conditions give the critical values of these coefficients, and 
substituting these critical values back into (43), one gets the stability boundaries. It is noted in [s] that in order 
to get results accurate to order g^, it suffices to work with Tq and Ti. Accordingly, we will ignore T2 and higher 
order "scales" and let x{t) = x{Tq,Ti). 



The equations. Using the notation Di = we have, up to second order in q, 

DqXo + q{DiXQ + Df^xi) + q^{DiXi + 0^X2) 



dx 
dr 

dT^ 



= Dlxo + q{Dlxi + 2D0D1X0) + q^{Dlxo + 2DqDiXi + 0^X2) ■ 



(44) 
(45) 



Substituting (44 1, (45) and (43) (where ag is arbirary, for now) in (38), we get, 

D^xo + q{Dlxi + 2DqDiXq) + q^{DlxQ + 2DqDiXi + D^xs) 

+ (ao + qai + q^a2)[xo + qxi + q^X2) + 2q{xn + qxi + q^X2) cos2To = . 

Collecting terms with similar powers of q and setting them to zero, we get the following partial differential 
equations: 



DIxq 



aoxa = 
agXi = —2D^DiXi^ 
0^x2 + aQX2 = -2DqDiXi 



Dlxi 



aiXQ — 2xq cos 2Tq 
Dlxo 



(46) 
(47) 
(48) 



0.QX2 = —'ZUqUiXi — U^xq — aixi — a2Xf) — 2xi cos2To 
As mentioned above, if the value of a under consideration is away from a set of critical values, one can ignore ai 



and 02, set a = ao, and solve the equations (46)- (48) to obtain an approximate solution to (38). This approach. 



in fact, is how the critical points are obtained in the first place: we set a ^ uq, and obtain the special values of 
flo for which the approach fails (due to the appearence of "small denominators"). See, e.g., [sjjo] for a discussion 
and examples. 

Here, we will skip this step, and just borrow the result that ao — and oq = 1 are the two critical values 
relevant for the first stability region near the origin. We next perform expansions around these two points. 



5.1.1 Expansion around a = 



Setting ao = in (46)- (48), we get 



D'oxo = 

DqXi — —2DoDiXo — aiXQ — 22;oCOs2To 

DqX2 — —2DqDiXi — D\xq — 02^0 — aiXi — 2xi cos2To . 



The general solution to ( 49 ) is 



xo{To,Ti) = A{Ti) + ToB{Ti) , 



(49) 
(50) 
(51) 

(52) 



where A and B are arbitrary functions. In order to suppress the secular term, i.e., the term that grows in time, 
we set _B = 0. This gives, 

a;o(To,Ti) = A(ri). (53) 



Substituting ( 53 1 into ( 50 ) , we get 



Llgxi = Ao(Ti)(-ai-2cos2ro). (54) 
The solution for xi has a secular term proportional to ai; in order to suppress unbounded growth, we set ai = 0. 



Solving ( 54 ) with this assumption and setting to zero another secular term that appears (this term is analogous 
to the B term in (52)), we get, 

xi{To,Ti) = C{Ti) + ^A(Ti)cos2To. (55) 
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Substituting (53 1 and (55) in (51) and using ai — 0, we get, 

Dlx2 ^2A'{Ti)sm2To~ A"{Ti) ~ a2A{Ti) - 2C{Ti) cos2To 

-^-^cos4To. 
2 2 

In order to avoid terms in X2 that grow linearly in Tq, we must have, 

A"(ri) + (a2 + ^)A(Ti) = 0. 

This equation will have non-growing solutions for A{Ti) if and only if, 

1 

Combining the conditions ao = 0, ai = 0, and a2 > ^1/2, we see that the curve in the a-q plane separating 
bounded solutions from unbounded ones around a = 0, g = is given by. 



1 2 



(56) 



with a > —q'^/2 being the stable region. 
5.1.2 Expansion around a = 1 



Setting ao = 1 in (|46|)-(|48|), we get, 
D^xo + xo 







DqXi + xi = — 2L)o£'ia;o ^ flia^o ~ cos 2ro 

0^X2 + X2 = —2DoDiXi — D\xq — aixi — a2a;o — 2x1 cos 2ro ■ 



The general solution of (57 1 is, 



Xo(To, Ti) = A(Ti) cosTo + B{T^) sinTo , 



Plugging this in (58), we get, 

Dlxi + xi =2A' sin To - 2B' cos To — ai{A cos To + -B sin To) 
- yl(cosTo + cos3To) + S(sinro - sinSTo) . 



(57) 
(58) 
(59) 



(60) 



(61) 



The solutions of this equation will have growing parts unless the coefficients of the "resonance terms" sin To 
and cos To on the right hand side vanish. This gives. 



2^' + (l-ai)B = 
2B' + {l + ai)A = 0. 



(62) 
(63) 



This system will have exponentially growing solutions if < 1, and oscillatory solutions if > 1. Thus, the 
critical values of ai are ai = ±1. Assuming > 1, the general solution is, 

A{T,) 



c sin ATi + d cos ATi 
2A 



fli — 1 



(-dsinATi +ccosATi) 



(64) 
(65) 



where A = y '-^ and c and d are constants. After substituting (|64|)-(|65| in (|6l|) and setting the resulting 
resonant terms to zero, we get, 

Dlxi +xi = -(Acos3To + Ssin3To), (66) 
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whose general solution is, 



a;i(ro,Ti) =C(ri)cosTo + i^(ri)sinTo 



A(Ti) 



cos 3Tn 



8 



— sinSTn. 



(67) 



Plugging (67) in (591, the right hand side of (591 ends up having various terms proportional to sinTp and cosTq- 
Once again, these resonant terms will result in growing solutions for X2, so we set their coefficients to zero. This 
gives. 



Recalling A" = -X^A and B'' 



2C' + (f-ai)D = B" + B{a2 + -) 

o 

2D' + {l + ai)C - -A" - A{a2 + \). 

o 

-- —X^B from above (see ( 641(65 1), we get, 
2C' + il~ai)D = S(i— ^+a2 + ^) 



2D' + {l + ai)C = -A{^—-^ 



0-2 



(68) 
(69) 

(70) 
(71) 



Combining these equations, using the explicit solution (64][65]) for A and B and once again requiring the absence 
of resonant terms, we get, 

a2 = -- + ^— . (72) 



Finally, combining the conditions for oq, ai and 02, we get the formula for the stability boundaries near q = 0, 
a = 1 as, 



a = ao + aiq + a2q 



(73) 
(74) 



with the region between the two curves corresponding to unstable solutions, and the region outside corresponding 
to stable ones. The results (56) and (74) are in agreement with the classical results for Mathieu equation. 
For future reference, we note that if the a in the Mathieu equation (38) was replaced with —aa, the stability 

(75) 



boundaries ( 56 ) and ( 74 ) would be replaced with. 



-1 
a 



1 , 

1±9 - 



and , 



1 

2a' 



(76) 



5.2 Two-variable, coupled Mathieu's equations 



We next turn to the coupled Mathieu system (20)-(21 ), which we reproduce here: 

i + ax + 2g(ca:; + sy) cos 2r = 
?/ — aay + 2g(sa; — cy) cos 2r = 0. 



(77) 
(78) 



These equations describe the radial motion of an ion in an asymmetric trap with decoupled axial motion. Recall 
that the coordinates are chosen to be along the DC principal axes, and the constants c and s are given in terms 
of the relative angle Q between the RF and DC axes as. 



c = cos 26* 
s — sin 2Q . 



(79) 
(80) 
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As in the single variable case, we begin the multi-scale analysis by promoting x and y to functions of two 
separate time variables, Tq — t and Ti — qr, and expanding them to series in q, 



X = xa{TQ,Ti) + qxi{Ta,Ti) + q^X2{To,Ti) 
y = yo{To,Ti) + qyi{T(i,Ti) + q^y2{TQ,Ti) 



If a is near a critical value ao where this expansion fails to work, we expand a as well, 

a = ao + aiq + a2q^ + . . . . 
Working to second order in q, and collecting terms with similar powers of q, we get. 



,2, 
0" 



D^xq + aoxo 



=0 



- a^xi = — 2DoDiXo — aixo — 2(cxo + syo) cos2ro 
>qX2 + aoa;2 = - 2DoDiXi - DiXq - aiXi 

— a2XQ — 2{cxi -\- syi) cos 2To 



for (77), and 



-DoJ/i ^ "«oyi = - ^DoDiyo + aaiyo - 2{sxo - cyo) cos2ro 
Dly2 - aaoy2 = - 2DoDiyi ~ Djyo + aaiyi 

+ aa2yo - 2{sxi - cyi) cos 2To 



(81) 
(82) 



(83) 

(84) 
(85) 

(86) 

(87) 
(88) 

(89) 



for (78). Following our discussion in Section 3.2 we will assume a > 0. For a > 0, the critical values are 
a — and a — 1 (once again we will skip the justification for these values, and refer the reader to ^). The 
expansion around a = will proceed similarly to the single variable case, but the a = 1 case requires special 



attention. Assuming oq > and a > 0, we see that the solutions of (84) are oscillatory, but the ones for (87) 
are exponential. 

Assuming that the initial conditions of the system are fine-tuned so that the unbounded solutions are not 
excited, one can investigate the bounded part of the general solution. However, as mentioned in our discussion 
in Section [3.2[ noise will undoubtedly excite the unstable solutions in practice. 

Nevertheless, our stability plots obtained by the techniques of Section [3?2| suggest that the curves separating 
regions of partial stability from those of full instability may in fact be of some use. We observe from the stability 
plots of Figures [2] and [s] (and from additional plots we will present below) that the primary stability region 
around the origin is bounded by an extension of such curves. In other words, the curves emanating at the 
critical point ao = 1 begin as boundaries between regions of full instability and regions of partial instability, 
but once they intersect with the stability boundaries emanating from ao = 0, they become boundaries between 
full instability and full stability. Motivated by this empirical observation, we will obtain formulas for the curves 
emanating at the critical point ao = 1, and use them as stability boundaries after they intersect with the curves 
emanating from ao = 0. We will demonstrate the reliability of this approach in the stability plots we present 
below. 



5.2.1 Expansion around a = 



We begin with the expansion around a = 0. Setting ao = in equations (84|-(89), we get, 

12„ 



DqXi = — 2Df)DiXo — aiXf) — 2cxo cos 2To — 2syQ cos 2To 



Dlx2 = - 2DaDiXi ~ Dfxo 
— 2syi cos 2ro , 



aiXi — a2Xo — 2cxi cos2ro 



(90) 
(91) 

(92) 
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and, 



^Ivi = - 2DoDiyo + aaiyo 
Dly2 = - 2DoDiyi - Djyo H 
— 2sxi cos 2To . 



- 2cyo cos 2Tq — 2sxo cos 2To 
aaiyi + aa2yo + 2cj/i cos 2Tq 



Solving (90) and (93 1 and setting the secular terms to zero, we get, 

xo{To,Ti) 



yo(To,Ti) 



A(Ti) 



Substituting these in (91) and (94 1 gives, 



D^xi = ~aiA{Ti)-2{cA{Ti)+sC{Ti))cos2To 



aaiC{Ti) + 2(cC(Ti) - sA(Ti)) cos2ro . 



(93) 
(94) 

(95) 



(96) 
(97) 



(98) 
(99) 



In order to avoid growing xi and j/i, we need to pick oi — 0. Solving (98) and (99|) with this assumption and 

^row linearly in Tq to zerc 

(cA(ri) + sC(ri)cos2ro 



setting the coefficients of two other secular terms that grow linearly in Tq to zero, we get 

1 



a;i(To,ri) 
yi{To,Ti) 



isA{Ti) - cC{Ti)) cos2Tq . 



(100) 
(101) 



Plugging the solutions (96 1, (97), ( |100 ), (101) in ( [92| and (95), we get equations for X2 and y2- These will have 
bounded solutions only if the To-independent terms on the right hand sides of the equations add up to zero. 
Asserting this condition gives. 



A"{T,) 



s') A{T,) 



_aa2 + -(c2 + s2))C(ri) 



(102) 
(103) 



Equations (102) and (103) will have oscillatory solutions for A{Ti) and B{Ti) when a2 > — l(c^ + s^), and 



simulatenous stability occurs for — 5 < ^2 < 
given by. 



In 



when a2 < ^{c^ + •8^)1 respectively. Using (79) and (80), we have, + = 1. Assuming a > 0, we see that 



other words, the stability boundaries around a = are 
-\q^ (104) 



Note that these conditions are what would be obtained from separate stability analyses of (77) and ( [78| , 
respectively, if the coupling terms in those equations were ignored and the single- variable results ( 56 1 and ( 76 ) 
of Sect ion fS . 1 . 1 1 were used. 

5.2.2 Expansion around a — 1 

1. As mentioned above, in this case we will seek solutions with partial stability. 



We next expand around a = 
For ao = 1, (p4l)-([86l) become. 



DIxq 
Dlxi 
Dlx2 



xo = 

xi — —2DoDiXq — flixo — 2{cxo + syo) cos 2To 
X2 = —2DqDiXi — D\xq — axX\ — a2Xo 
— 2{cxi + syi) cos2To , 



(106) 
(107) 
(108) 
(109) 
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and (87l-(89) become, 



D^yo - aya =0 
D^yi -~ayi=- 
Dly2 - a?/2 = - 



2DoDiyo + aaiyo - 2(sa;o - cyo) cos 2To 
2DoDiyi - Dfyo + aa2yo 

+ aaiyi — 2{sxi — cyi) cos 2To . 



The general solution to ( 106 1 is, 



a;o(To, Ti) = A(Ti) cos To + sinTo , 



Assuming a > 0, the solution to (110) is exponential, 



2/o(To, Ti) = S(Ti) exp (V^Tq) + F(Ti) exp (-v^Tq) . 



(110) 
(111) 

(112) 



(113) 



(114) 



This shows that the general solution to the coupled system is unstable for a > 0, cq = 1, as discussed above. 
In order to investigate partial stability, we set the coefficients of the exponential solutions to zero, E{Ti) = = 
F{Ti), which gives ?/o(2o,Ti) — 0. Substituting this and (1131 in (107), and setting the coefficients of sinTg 
and cos Tq on the right hand side to zero in order to avoid resonant (secular) terms that would result in growing 
solutions for xi, we get. 



Combining these, we get. 



2A' - (ai - c)B 
2B' - (ai + c)A 



A'' 



(115) 
(116) 



-A 



-B. 



(117) 
(118) 

which will have oscillatory solutions if of > (? and exponential ones if a\ < (? . Thus, the critical values for 
partial stability are, 

ai=±c. (119) 

The oscillating solutions for A and B are given as. 



^m) 

i3(Ti) 



r sin \T\ + p cos \T\ 
2\ 



fli — c 



(-p sin ATi + r cos ATi ) , 



(120) 
(121) 



where A — 



-, and r and p are constants. 



After having ensured that the resonant terms in (107) vanish by enforcing (1151 and (1161, we can solve 



(1071 to get the oscillatory solutions for x\. We get, 

a^i (To , Ti ) =C(Ti ) cos T^^D (Ti ) sin Tq 

+ ^A(Ti) cos3To + ^B(Ti) sin3ro . 



(122) 



Similarly, substituting the solution (113) in (111), we get a bounded solution for y^ after setting the coefficients 
of the exponential terms to zero: 



2/i(7o,Ti) 



(a + 1) 



(yl(ri)cosTo-B(ri)sinro) 



(a + 9) 



(A(ri) cos 3ro + B(Ti) sin 3ro) 



(123) 
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We next substitute (122) and (123) into (108) and collect the resonance terms, i.e., terms proportional to sinTo 
and cos To, on the right hand side. Setting the coefficients of these resonances to zero in order to avoid growing 
solutions for X2, we get. 



where. 



Using (117)-(118l, we get 



2C'(Ti) + i?(ri)(-ai+c) = B{Ti)" + I3B{T^) 
2D\T,) + C{T,){ai+c) = -A''{T,)~(3A{T,). 



where ^ = /3 — = a2 + 
the decoupled equations, 



0-2 



c2 2s2(5 + a) 
8" ^ (9 + a)(l + a) 



2C'{Ti) + D{Ti){-ai+c) = p.B 
2D'{T,) + C{Ti){ai+c) = -/iA , 



(124) 
(125) 

(126) 



(127) 
(128) 



2s^(5+a) 
(9+a)(l+Q) 



. Combining (127) and (128) with (1151 and (116), we get 



C" 
D" 



Qi - = 

4 2 



(129) 
(130) 



Now, since the solutions (1201 and (121) for A and B are in resonance with the left hand sides of (129) and 



(130), in order to avoid growing solutions. 



the coefficients on the right hand sides of these equations must vanish. Using (119), this gives. 



02 = 



8 (l + a)(9 + a)' 

Thus, the second order approximation to the relevant stability boundary starting at a = 1 is given as. 



1 — cq ■ 



c2 2s2(5 + a) 
¥ ^ (l + a)(9 + a) 



(131) 



(132) 



where we picked the negative sign for the first order term in order to get the curve that approximates the upper 
boundary of the primary stability region. 

Boundary of the primary stability region. In order to obtain the approximate boundaries of the primary 
stability region, we also need to find an approximate formula for the stability boundary that emanates from 
negative a when q = 0. Our expansion around a = 1 was based on the assumption that a > 0, so that x had 
oscillatory solutions to zeroth order, and y had exponential ones. In order to obtain the curves for negative a, 
we just replace the roles of x and y by a series of redefinitions. Namely, we set. 



x^y, y^x, a 



-aa , a = l/a, c. = —c, s 



(133) 



These redefinitions transform the equations of motion ( 77 ) and ( 78 ) into exactly the same form, with the 



variables (except q and t) being replaced by their tilded counterparts, 

X + ax + 2q{cx + sy) cos 2r = 
y ~~ aay + 2q{sx — cy) cos 2t = 0. 



(134) 
(135) 



Now, if a < 0, the transformation (|133|) makes a > 0, thus, we can apply (132) to (134)-(135) by replacing all 

a^l-cq-[— + 7^^,^ ) q' , (136) 



the quantities in ( 132 ) by their tilded counterparts. This gives, 

2S2(5 + a) 



8 (l + a)(9 + (5) 
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or, equivalently, 



[l-cq- 

a 



8 



2s^{5 + l/a) 
(l + l/a)(9 + l/a) 



(137) 



This gives the required approximate boundary of the primary stabihty region for a < 0. Note that the curve 
starts at a = — l/a. 

In Figures |12|18t we compare the results of this section with the results from the numerical analysis of 
Section 3.2 for various values of a and 



6 Conclusions and discussion 



In Figures[T2] THj we present the results of our stability analysis. We compare the results of the numerical method 
(purple areas) to the stability boundaries obtained by the multi scale perturbation analysis (black 



of Section 3.2 



lines). We show plots for a range of angles 6 between the RF and DC axes, and a range of as (recall equations 



(20l-(21|). A few important conclusions are evident from the plots. 



• The primary stability region does not get smaller when a nonzero angle 6 is introduced between the RF and 
DC axes if the other variables are kept fixed. Such a nonzero angle between RF and DC fields represents 
the case of an asymetric surface electrode geometry and/or asymmetric voltages. 

• Although the primary stability region does not change appreciably when 9 is varied, two secondary stability 
regions are highly variable, and when 9 has the special value of 45°, they join the primary region to result 
in an exceptionally large region of stability. 

• The curves obtained from multi-scale perturbation theory approximate the boundaries of partial stabilit}!^ 
quite accurately when q is near zero. Regions of partial stability become regions of full stability when they 
"overlap" near the primary stability region of classical, symmetric traps. Unfortunately, in this region, q 
is likely large enough to make the accuracy of the approximation unsatisfactory, especially for angles 9 
close to the critical value of 45° . Since, in practice, partial stability means instability, the small q region, 
where our multi-scale formulas are uniformly accurate, is not of relevance to the practical problem of the 
stability of ion motion. 

• When we lay the approximate stability boundaries for the symmetric Paul trap (described by decoupled 
equations) on top of the coupled stability plots, we see that these boundaries consistently underestimate 
the size of the stability region for the coupled system. 

A relative angle between the RF and DC axes does not change the boundaries of the primary stability region 
significantly, unless the angle is close to 45°, in which case the stability region is enhanced. However, the shapes 
and locations of a pair of secondary stability regions are highly dependent on the relative angle. 

The conclusion of this analysis for the practical problem of asymmetric ion trap design and operation is as 
follows: Just as in the case of symmetric Paul traps, the q-a stability plot of the standard single- variable Mathieu 
equation is sufficient to allow one to determine stable operating conditions for asymmetric surface traps with 
long RF electrodes. It is possible to proceed as in the case of symmetric Paul traps, by obtaining a pair of "a" 
values for the two radial principal axes of the DC potential, and a pair of "g" values for the two radial principal 
axes of the RF potential. Ignoring the fact that there is a nonzero relative angle between the RF and DC axes, 
the two decoupled Mathieu equations for these {q, a) pairs can be used to determine the stability properties of 
the coupled system. This approach does not give precise stability boundaries for asymmetric traps, but it is 
"safe" in the sense that trap operating conditions deemed stable by this method will in fact be stable for the 
coupled system. By ignoring the coupling, the size of the primary stability region is simply underestimated. 



^'^In Figure [12] we do not show the plot for 9 = 90° since it is the same as the plot for 9 = 0°. In fact, as can be seen from the 
figure, the stability behavior of the system is symmetric around 9 = 45° , in the sense that the plots for 45° -I- A9 and 45° — A9 
are identical. This is due to the fact that the equations of motion are invariant under the transformation 9 — > it/2 — 9, y —y. 
Since the system is symmetric under reflections of the y axis, this transformation should leave the stability plot invariant. In 
Figures |l3|18[ we only show the stability plots for 9 between 0° and 45°. 

'^^If there is a set of initial conditions for which the ion motion is bounded and another set for which it is unbounded, we call the 
system partially stable. Since it is impossible to tune the system perfectly, the unbounded solutions of a partially stable operating 
point of the trap will get excited in practice, and the ions will be lost. 
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Figure 12: Stability plots for a = 0.5 and 9 between 0° and 90°, with a step size of 5.625°. The x-axis is q and 
the y-axis is a in each plot. The areas shown in dark purple correspond to complete stability, and those in light 



purple correspond to partial stability, predicted using the numerical method of Section 3.2 The black curves 
are approximate stability boundaries for the coupled, two variable Mathieu system, obtained from equations 



(104 1, (105), (132), and (137). The red curves are the approximate stability boundaries for the corresponding 



decoupled systems, obtained from equations (74) and (75). The approximate boundaries around a — for 



the coupled and the decoupled systems are identical (see equations (56)-(76) and (104)), thus, we only show 
black curves for the approximate boundaries passing through a = 0. (see the final paragraph of Section [6] for a 
discussion of the relation between the coupled system and the decoupled system) . We see that the black curves 
follow the partial stability boundaries quite accurately when q is small, however, when 9 gets close to 45°, the 
curves starting at a = 1 and a = —l/a = —2 (the latter point is outside the region shown in the plots) lose their 
accuracy near the primary region of full stability. The colored curves representing the approximate boundaries 
of the decoupled system, while inaccurate predictors for the coupled system, consistently underestimate the size 
of the primary stability region. 
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Figure 18: Stability plots for a — 2.5 and 9 between 0° and 45°, with a step size of 5.625°. See the discussion 
under Figure 12 
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